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ABSTRACT 

The Energy Finite Element Analysis (EFEA) has been utilized successfully for modeling complex structural-acoustic systems with 
isotropic structural material properties. In this paper, a formulation for modeling structures made out of composite materials is 
presented. An approach based on spectral finite element analysis is utilized first for developing the equivalent material properties for 
the composite material. These equivalent properties are employed in the EFEA governing differential equations for representing the 
composite materials and deriving the element level matrices. The power transmission characteristics at connections between members 
made out of non-isotropic composite material are considered for deriving suitable power transmission coefficients at junctions of 
interconnected members. These coefficients are utilized for computing the joint matrix that is needed to assemble the global system of 
EFEA equations. The global system of EFEA equations is solved numerically and the vibration levels within the entire system can be 
computed. The new EFEA formulation for modeling composite laminate structures is validated through comparison to test data 
collected from a representative composite aircraft fuselage that is made out of a composite outer shell and composite frames and 
stiffeners. NASA Langley constructed the composite cylinder and conducted the test measurements utilized in this work. 

INTRODUCTION 

The EFEA is an emerging technology for analyzing complex structural-acoustic systems [1-18]. The EFEA utilizes the spatially 
averaged vibrational or acoustic energy as a primary variable for developing governing differential equations for each wave type 
present within a system. The fundamental solution to the governing differential equations presents an exponential variation with 
space, which eliminates the requirement of using small finite elements at high frequencies. At discontinuities (i.e. connections 
between structural components, structural and acoustic interfaces, changes in material properties, changes in geometry, etc.) the 
primary variables in the EFEA become discontinuous. Joint matrices must be developed in order to introduce the power transfer 
characteristics at each discontinuity into the EFEA system of equations. For each finite element encountered at a discontinuity, the 
joint matrices represent relationships between the power flow at the element’s boundaries and the corresponding EFEA primary 
variables. The joint matrices are developed from power transfer coefficients representative of the connections between the members 
encountered at the discontinuity. The derivation of power transfer coefficients is based on considering a connection between semi- 
infinite members, prescribing an impinging wave type originating from one of the members, and evaluating all transmitted and 
reflected wave types in all members. A library of element matrices and a library of joints can be combined for creating the EFEA 
model for a complex system and computing its structural-acoustic response. 

This paper aims to further extend the capability of the EFEA method for computing the vibro-acoustic response of aircraft type of 
structures made out of composite laminates and/or composite sandwich panels. A Spectral Finite Element Method (SFEM) [19, 20] is 
used for evaluating the group speed Cg and the damping rj* due to the effectiveness of SFEM for taking into account the layer- wise 
composition of composite panels along with the shear deformation effects. The c* and r]* values computed by the SFEM are used in 
the EFEA formulation for deriving the element level matrices of the composite materials. A method for deriving the power 
transmission coefficients between composite members is presented. The new method can take into account the shear deformation 
effects of composite plates. The power transfer coefficients are used for deriving the joint matrices in the EFEA formulation. 
Theoretical aspects of these two formulations are presented in this paper. A numerical EFEA implementation for composite materials 
is utilized for analyzing the NASA Composite Test-bed Cylinder (CTC) and comparing simulation results to test data for the vibration 
at several bays of the cylinder. The CTC has longitudinal stiffeners (representing stringers) and circumferential stiffeners 
(representing frames) and represents the dynamic characteristics of a composite aircraft fuselage [22]. Excitation is applied through a 
shaker at different locations on the CTC and measurements are collected using a laser vibrometer at twelve bays of the cylinder. The 
results from the EFEA simulations are compared to test data in the frequency range that covers the 1/3 octave bands of 800Hz - 
4,000Hz. The experimental set-up and the correlation study are discussed in this paper. 
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DERIVATION OF GROUP SPEED AND DAMPING FOR COMPOSITES 


The Cg and 77 * values computed by the SFEM are used in the EFEA formulation for composites for deriving the element level 
matrices. The basic steps of the SFEM approach are presented here. Figure depicts a plane wave propagating in the positive x 
direction with frequency co and wavenumber k in a multilayer composite plate. The x direction comprises a reference for orienting the 
various layers. The through-thickness discretization in the z-axis along with the assumption of the harmonic wave motion in the % 
direction results in the following form of displacement field, u r = [u, v, w] at any point, (x, y, z ) within the plate, 


u(x,z,t) = N^ue 1 ^ kx>} ( 1 ) 

where N (z) is a matrix of shape functions and u is a vector of the nodal degrees of freedom of the form: 

u =[u 1 ,v 1 ,w 1 u 2 , v 2 , w 2 ••• u Ne+1 ,v Ne+1 ,w Ne+1 ]T (2) 

where N e is the total number of linear finite elements through the thickness. It is noted that any of the displacement components in u 
may be complex numbers. 



Figure 1 : A plane wave propagation in a multilayer composite panel 
The time-averaged total kinetic and potential energies, (T) and (V) are given by 

(T) = ^-/ n pu H udQ.; (7) = En ( 3 ) 

where H stands for the Hermitian transpose. Substitution of appropriate stress-strain and strain-displacement relations and the 
displacement field equation (1) into equation (3) results in: 

(T) = u h cl) 2 Mu; ( V ) = u H Ku (4) 

where K and M are, respectively, the stiffness and mass matrices. The replacement of the spatial derivatives with respect to x and z by 
- i k and dN/ dz can yield the expressions for the stiffness and mass matrices. Using Hamilton’s principle results in: 

[K(fc) - co 2 M]u = 0 (5) 

Since K = K(/e), for a prescribed circular frequency co, equation (5) can be solved for the eigenvalues co 2 and eigenvectors u. The 
hysteretic damping model can be applied to each layer to yield the following form for the total time-averaged energy loss associated 
with an arbitrary wave type, 

N e 

(Ediss) = 'Y J VMe)i ( 6 ) 

1 = 1 

where r]i and {e) t are, respectively, the structural loss factor and energy density of the Ith layer of multi-layered composites. Since 
(e) = (T) + (V) = 2(F), the damping loss factor associated with each propagating wave can be expressed as 
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( 7 ) 


Y.i= 1 V l u w K®u 

V= u tf Ku 

Here, the distinction should be made that K® is the stiffness matrix of the ith layer and K is the assembled global stiffness matrix of 
all layers. Equation (7) is derived for each wave heading angle and thus the angle-average must be considered for accounting for 
the full range of wave propagation directions when deriving the angle averaged damping loss factor as follows: 


„ fg V (.0i)dGi 
T] = — 

fg de i 


( 8 ) 


The variation of the angle 6 t is obtained by rotating the relative orientation of the x-axis with respect to the composite material (Figure 
1). The circular frequency for the wave number k can be written in terms of the Rayleigh quotient as 

u w K(fc)u 

" = TSmT (9) 

Figure 2 shows the wavenumber of an arbitrary wave propagating in an anisotropic medium as a function of wave heading, 0 t . The 
distribution and the direction of energy flow in anisotropic media like composites exhibit directional dependence, thus the direction of 
the energy flow, depicted as 0 e in Figure 2, is usually different from the wave heading, 0 t . Thus, an appropriate correction must be 
introduced in the calculation of the group speed. 



The group speed in the direction of wave propagation is c gd = dco/dk , and the differentiation of equation (9) with respect to the 
wavenumber, /c, results in: 


_ u H (dK(k)/dk)u 
Cge ~ 2a)u H Mu 


Referring to Figure 2, the following relationship holds: 


c _ c 9 e 
9 cos (Qi - Q e ) 


( 10 ) 


( 11 ) 


where the heading of group speed, 6 e may be derived from the geometric interpretation of the wavenumber curve shown in Figure 2 


tan 6 P = — _ 

dky/dOi 


dk x /d0i (d/c(0i)/d0i)cos#i - k^d^smOi 


v y/ wi (d/t(0i)/d0i)sin0j + /c(0j)cos0j 

Therefore, the angle-averaged group speed can be evaluated for the full range of wave propagation directions as follows: 


( 12 ) 
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J e c fl(0j) d0 i 

J 0 de t 


( 13 ) 


In the EFEA formulation for composites, c* g and rj* from Equations (13) and (8) respectively, are used when deriving the element level 
EFEA matrices. 


POWER TRANSFER COEFFICIENTS BETWEEN COMPOSITE MEMBERS 

The derivation of power transfer coefficients is based on considering a connection between semi-infinite members, prescribing an 
impinging wave type originating from one of the members, and evaluating all transmitted or reflected wave types in all members. 

Once available, the power transfer coefficients are used for computing the EFEA joint matrices which connect elements across a 
discontinuity in the EFEA system of equations. The formulation for computing power transfer coefficients between composite 
members is presented here. It must be noted, that this formulation accounts for First order Shear Deformation Theory (FSDT) effects, 
which are important when considering sandwich composite panels. The extended wave dynamic stiffness matrix approach [21] is used 
to calculate the wave power transfer coefficients of coupled composite members. 

Figure 3 presents a junction consisting of N semi-infinite composite members. The coordinate systems and the translational and 
rotational displacements are also presented in Figure 3. 



(a) 


z r w 



Figure 3: (a) A general N-plate junction and global coordinate system; (b) local coordinate system and displacements for plate n 

For either a laminated or a sandwich composite member, the FSDT allows writing the dynamic equations of motion in the following 
matrix- vector form: 
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Lu = 0 


(14) 


where u = [u, v, w, cp x , cp y ] T and L is a linear differential operator which can be expressed in terms of elastic constants A tj - , B and Z^y 
and inertial properties I a (a = 0,1,2). Substitution of a wave form, u = \\e~ lkxX+?iy+lcot into equation (14) yields 

(A 2 b 2 + Abi + b 0 )u = 0 (15) 

where b 0 , b 1? and b 2 can be obtained from L by replacing ^ and by (-i/c x ),A, and co. This quadratic eigenvalue problem can be 
solved for five pairs of A’s and u’s for given k x and a). Then, the edge displacements, d y and the corresponding elastic tractions, t y of 
each plate can be expressed in terms of the reflected/transmitted wave amplitudes c = [c lf c 2 , c 3 , c 4 , c 5 ] T and an incident wave 
amplitude c 0 as follows 


dy — [dyl dy 2 

dy5 ] C + d y0 Co 

(16) 

ty - [tyi *yz 

ty 5 ]c + ty 0 C 0 

(17) 


where dy^and t y j(/=l,2,. . .,5) represent edge displacements and elastic tractions due to the zth transmitted or reflected wave and d y0 
and t y0 are those for incident waves. Equations (16) and (17) can be combined to give 


ty = Kdy-f y0 (18) 

where K = [d yl d y2 ••• d y5 ] _1 [t yl t y2 ••• t y5 ] is the wave dynamic stiffness matrix and f y0 = Kd y0 - t y0 is the force due to the 
incident wave. It should be noted that f y0 = 0 in the plates where no incident wave exists. By implying subscript “y” for the 
remaining discussion, equation (18) can be rewritten for the nth plate as 


K n d n 


(19) 


The wave dynamic matrix is assembled to yield an equation for the computation of power transmission coefficients. Since d n and t n 
are defined in the local coordinate system, they need to be transformed to those at the common junction defined in the global 
coordinate system, dy and ty as follows: 

d n = R n d ; ; t j = R T n t n (20) 

where R n is a simple coordinate transformation matrix consisting of cosi p n and sim p n . Since the force equilibrium equation for the N 
plates connected through a structural joint is £n=i R^t n — 0, using equations (19) and (20) results in: 

(21) 


r N i 


r N l 

y R?XR n 

d J = 

Y Rjfn 

-n=l 


-n=l 


where dy can be evaluated from Equation (21). Using equation (16) the wave amplitudes are evaluated: 

C = [dyl dy 2 ••• dy 5 ] 1 (dy — dy 0 C 0 ) 


Then, the power transfer coefficients can be readily calculated from: 


r 0 -(o),60 = 


2 


I c 0i I 


2 


( 22 ) 


(23) 


where Cy is the amplitude of a reflected or transmitted wave type j for the given c oi , which is the amplitude of an incident wave type i 
with frequency co and an incidence angle of 6 . 


The above computed transmission coefficients are highly dependent on the angle of incidence, 6 for anisotropic media (i.e. composite 
panels). Thus, the power transfer coefficients are obtained through averaging of the power transmission coefficients with respect to the 
angle of incidence: 
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(24) 



T ij(fi) c gy(8) / 

c(e)c ge (e) “7 , 


Cgyje) 

c(O)c ge (0) 


dB 


The power transfer coefficients are used for creating the EFEA joint matrices which connect elements across a discontinuity. 
Composite material properties can be considered for the components connected across a discontinuity based on the formulation 
presented in this Section. 


SIMULATIONS AND COMPARISON TO TEST DATA 


The EFEA capability for modeling the response of composite aircraft type structures in high frequencies due to structure-borne 
excitation is validated by comparing simulation results to test data collected using the NASA CTC. Test results are compared to 
EFEA simulations in the frequency range that covers 1/3 octave bands from 800Hz - 4,000Hz. The description of the test set-up, and 
the correlation between the EFEA and the test data are discussed in this Section. 



(a) 



(b) 

Figure 4: Exterior (a) and interior (b) of the NASA composite test-bed cylinder 

The NASA CTC structure, presented in Figure 4, is a floor-equipped stiffened composite cylinder that is 1.676 m in diameter and 
3.658 m long. The 1.7 mm thick shell is made of carbon fiber filaments embedded in an epoxy resin. Specifically, the shell is made of 
9 layers of 0.189 mm thick T500-tape with a ply sequence of ±45, ±32, 90, +32, +45. The shell is reinforced by orthogonal supports, 
referred to as ring frames and stringers that are riveted and bonded to the skin. Ten J-section ring frames and twenty-two hat section 
stringers segment the cylinder into 198 rectangular bays with nominal dimensions of 0.203 m by 0.381 m. The entire structure is 
supported by two 88.9 mm thick particle board end caps. Details regarding the construction of the cylinder can be found in [22]. 

The structure was excited by a single electromagnetic shaker suspended from an adjustable excitation stand. The shaker was attached 
at ten different locations as indicated with red dots in Figure 5. Two sets of four excitation locations placed on the shell comprise 
groups A and B while excitation locations C and D are on axial stiffeners. The multiple excitation locations within groups A and B 
enabled spatial averaging to smooth out local variations in the driven bay. Multiple experiments were also conducted with the shaker 
attached at the same point to evaluate response sensitivity to unintentional variations in shaker alignment and preload. An impedance 
head was placed between the stinger and the structure to measure the axial acceleration and force at the drive point. The normal shell 
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velocity was measured using a laser vibrometer. The scanning laser vibrometer measured sidewall velocity at 204 points within the 
measurement area (shaded area) shown in Figure 5. 


Group B 



Figure 5: Drive point locations and measurement area (shaded area) 


The measured vibro-acoustic response of the test structure is compared to numerical results obtained using the EFEA formulation 
presented in this paper. The EFEA model of the NASA CTC is shown in Figure 6. The comparison between measurements and 
predictions is performed on the twelve shaded bays shown in Figure 5. Specifically, the time and space averaged energy density in 
each of the twelve bays is calculated from measurements and compared to the EFEA results. In order to facilitate the comparison, 
both the predictions and the measurements are normalized with respect to the input power, n iu . The final comparison is performed 
using the averaged normalized energy density, calculated as 


cLB = 10 log 


fy N C ( PW 2 T 


N r 






(25) 


where N c is the number of test cases (defined as the product between the number of repeated experiments and the number of excitation 
locations), (v) and (n in ) are the averaged normal shell velocity and the averaged input power. It should be noted that the energy 
density is related to the velocity squared as ( e ) = p(v) 2 . 



Figure 6. EFEA model for NASA CTC (without the end caps included in the figure) 
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Figure 7 shows the difference (in dB) between measurements and predictions in the 2000 Hz 1/3 -octave band. The twelve rectangles 
in each subplot correspond to the twelve measurement bays depicted in Figure 5. The four subplots correspond to the different 
excitation groups. In this case, positive values indicate that the measurements are higher than the predictions. In most cases, good 
correlation exists between the measured data and predictions with maximum discrepancies of 6 dB in select bays. Although not 
shown, similar results were observed in the other 1/3 -octave bands from 800 through 4,000 Hz. 



dB difference between EFEA and test results (Drive dB difference between EFEA and test results (Drive 
Point C) Point D) 

Figure 7: Difference (dB) in normalized energy density in the 2000 Hz 1/3 octave band 


Figure 8 presents the differences in the energy density between test and simulations for each panel, averaged over the entire frequency 
range between the 800Hz and 4,000Hz 1/3 octave bands for the four excitation locations. The frequency averaged difference between 
EFEA and test is computed as: 

Tijl^dBEFEA ~ dB test | ^ 

where, Nf is the number of center frequencies within the range of 1/3 octave bands from 800-4000Hz. The calculation of dB EFEA and 
dB test is performed according to Equation (25). The majority of the differences are below 3dB, which comprises a very good 
correlation. Further, in the panels where isolated 5dB differences are observed, the adjacent panels exhibit much better agreement, 
which indicates a lack of any physics based modeling errors. 
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dB difference between EFEA and test results (Drive 
Point Group A) 


dB difference between EFEA and test results (Drive 
Point Group B) 




dB difference between EFEA and test results (Drive dB difference between EFEA and test results (Drive 
Point C) Point D) 

Figure 8: Difference of normalized energy densities over measurement bays (in dB) 


SUMMARY 

A comprehensive formulation is presented in this paper for evaluating material properties and power transfer coefficients that can be 
used in an EFEA formulation for modeling composite structural systems. The new developments are employed for modeling a 
composite stiffened cylinder representative of a composite aircraft fuselage. Simulation results are compared successfully to test data, 
demonstrating how the new EFEA developments can be used for modeling a complex structural system made out of composite 
materials. 
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